#!/usr/bin/perl -w
=head USAGE

  panda_parse_HMMOUT_into_db.pl

=head PURPOSE

    This script takes the output files from HMMSearch and enters them
        into the specified postgres panda database. The original
        sequences used as the input into HMMSearch are also required

=head REMARKS

=cut
#die and warn check next

use strict;
use Getopt::Long;
use warnings;
use DBI;
use DBD::Pg;
use panda_parse_fasta;
use String::CRC;
use POSIX;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    panda_parse_HMMOUT_into_db.pl
                            --seq_file        <sequence.fasta>
                            [--dbname        <panda>            ]
                            [--user           <any_db_username>  ]
                            [--password       <secret>           ]
                            [--help]
							<HMMFiles ...>

  				
  				
    This script takes the output files from HMMSearch and enters them
        into the specified postgres panda database.
    The original sequences used as the input into HMMSearch are
        also required (<sequence.fasta>).

	
USAGE



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# options variables
# options variables
my $help;
my $dbname			='panda';
my $username		='';
my $userpassword	='';
my $seq_file		='';


GetOptions (
			'help'				=> \$help,
			'dbname=s'			=> \$dbname,
			'user=s'			=> \$username,
			'password=s'		=> \$userpassword,
			'seq_file=s'		=> \$seq_file,
			);



die $usage if ($help);
die "\nERROR:\n\tMissing Argument:(dbname)\n" .$usage
				unless ($dbname);
die "\nERROR:\n\tMissing Argument:(seq_file)\n" .$usage
				unless ($seq_file);
die "\nERROR:\n\tMissing HMM files:\n" .$usage
				unless (@ARGV);






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Various functions

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
sub parse_HMM_output_file($$\%$$\@\@);
sub save_domain_data_in_db($\@\%$);
sub get_dom_code_thresholds($$$$);
sub fasta_to_pi_callback($$$$);
sub drop_indices($$);

my $cnt_dom_prot;
#_________________________________________________________________________________________

#	drop_indices($$)

#		drops all indices for a given table and
#		returns the commands needed to recreate them

#_________________________________________________________________________________________
sub drop_indices($$)
{
	my ($dbh, $table_name) = @_;

	# get all index names and the commands to recreate them
	my $sth_get_indices	= $dbh->prepare("SELECT indexname, indexdef FROM pg_indexes "
							."WHERE tablename = ? ");
	$sth_get_indices->execute($table_name) or die "\nError:\n$DBI::errstr\n";
	my @recreate_commands;
	my @index_names;
	while (my @res = $sth_get_indices->fetchrow_array())
	{
		push(@recreate_commands, $res[1]);
		push(@index_names, $res[0]);
	}
	
	# return nothing if no indices on this table
	return () if (@index_names == 0);
	
	# drop all indices
	for (@index_names)
	{
		$dbh->do("DROP INDEX $_")
				  or die "\nError:\n$DBI::errstr\n";
	}
	
	return @recreate_commands;
}

#_________________________________________________________________________________________

#	parse_HMM_output_file

#		returns the number of HMMs which were used to produce this file
#_________________________________________________________________________________________
sub parse_HMM_output_file($$\%$$\@\@)
{
	my ($HMMoutput_file_with_path,
		$seq_file, $acc_to_pi, $dbh, $seq_cnt, $obsolete_hmms,
		$HMMs_sans_hits) = @_;

	# open file
	open HMMOUT, $HMMoutput_file_with_path or
				die "\nError!\nCould not open HMM output file [ $HMMoutput_file_with_path ]\n$!\n";


	
	$HMMoutput_file_with_path =~ m{.*/(.*)};
	my $HMMoutput_file = $1;
	my %dom_codes;

	# loop through each record in HMM
	while (!eof(HMMOUT))
	{
		#hmmoutput files begin with hmmsearch
		while (<HMMOUT>)
		{
			chomp;
			next if $_ eq '';
			unless (/hmmsearch/ || /hmmpfam/)
			{
				print STDERR "\n!!!! Major Error !!!!\n "
					."The hmmoutput file $HMMoutput_file record beginning on line $. ".
					"does not contain either 'hmmsearch' or 'hmmpfam'!!\nline $.:[$_]\n\n";
				# skip to end of record
				while (<HMMOUT>)
				{
					last if m{^//};
				}
				next;
			}
			last;
		}
		last if (eof(HMMOUT));

		# skip to and get HMM file name
		#	eg. "HMM file:  /share/databases/PFAM/123.hmm [???]"
		#                                                  ^^^^^^^
		my $HMM_filename;
		while (<HMMOUT>)
		{
			next unless /^HMM file: +(.*)$/;
			$HMM_filename = $1;

			# remove path, square brackets and spaces
			$HMM_filename =~ s{.*/}{};
			$HMM_filename =~ s{ *\[.*\] *}{};
			
			last;
		}
		
		# skip to and get HMM name and accession eg. "14-3-3|PF00244|14-3-3 proteins"
		#                                             ^^^^^^ ^^^^^^^
		my ($HMM_name, $HMM_acc);
		while (<HMMOUT>)
		{
			next unless /^Query HMM: ([^|]+)\|([^|]*)\|/;
			($HMM_name, $HMM_acc) = ($1, $2);
			last;
		}

		# skip one past the line with the word "Sequence"
		while (<HMMOUT>) {last if /^Sequence/}
		<HMMOUT>;


		# save a list of protein E-values
		my %protein_e_values;
		while (<HMMOUT>)
		{
			
			# IPI:IPI00000816.1|				  14-3-3 protein    527.1   1.4e-154   1
			# gi|479627|pir||S34968               fibulin, splice   108.3    1.9e-27   3
			last unless (/^([^:|]+[|:][^| ]+)	# gi or pi or IPI or swissprot
						.*\s+					# rest of accession, Description, Score
						([0-9eE\-.]+)			# E-value
						\s+\d+$/x);			# N
			$protein_e_values{$1} = $2;
		}


		# skip one past the line with the word "Sequence"
		while (<HMMOUT>) {last if /^Sequence/}
		<HMMOUT>;

		# list of proteins for this domain
		my %domain_data_per_pi;

		while (<HMMOUT>)
		{
			# Sequence                   Domain  seq-f seq-t    hmm-f hmm-t      score  E-value

			# gi|5803225|ref|NP_006752.1|   1/1  4   239 ..   1   236 []   525.1 6.3e-153
			# IPI:IPI00000816.1|            1/2  5   539 ..   1   126 []   310.1 9.3e-153
			last unless /^([^:|]+[|:][^| ]+)	# gi or pi or IPI or swissprot  !save! $1
							\S*\s+				# the rest of the accession
							\S+\s+				# the domain
							(\d+)\s+			# start  of the sequence		!save! $2
							(\d+)				# finish of the sequence		!save! $3
							/x;					#

			my ($acc, $b, $e) = ($1, $2, $3);
			my $pi = $acc_to_pi->{$acc};

			# check that pi could be retreived
			die "\nFatal Error!\n "
				."The protein sequence for\n\t[ $acc ] matched by the HMM\n\t[ $HMM_name ] "
				."was not found within\n\t[ $seq_file ]. "
				."Please supply the original sequence file passed to HMMSearch to "
				."generate\n\t[ $HMMoutput_file ].\n"  unless defined $pi;

			# check that the protein e_value exists
			die "\nFatal HMMOutput Parsing Error!\n".
				"\tThe per protein E-value for [ $acc ] on line $. of [ $HMMoutput_file ] could not be found.\n"
					unless exists $protein_e_values{$acc};
			

			# save domain data linked to pi
			$domain_data_per_pi{$pi} = [$b, $e, $protein_e_values{$acc} / $seq_cnt];
		}

		
		my $hmm_seqcnt;
		# skip to sequence count
		while (<HMMOUT>)
		{
#			last if m{^//};
			if (/Total sequences searched: +(\d+)/)
			{
				$hmm_seqcnt = $1;
				last;
			}
		}

		# make sure there is a sequence count
		unless (defined $hmm_seqcnt)
		{
			print STDERR	"\nWarning!: No count of searched sequences "
							."for $HMM_name [ $HMM_acc ] in [ $HMMoutput_file ]\n";
			next;
		}
		# make sure the sequence count is correct
		unless ($hmm_seqcnt == $seq_cnt)
		{
			die "\nError! The sequence count for the hmm output file [ "
				."$HMMoutput_file:$hmm_seqcnt ] is different from that for "
				."the supplied sequence file [ $seq_file:$seq_cnt ]\n";
		}

		# die if end of file prematurely encountered
		if (eof(HMMOUT))
		{
			print STDERR "\nError:\n".
				"Premature end after $. lines in HMM output file [ $HMMoutput_file_with_path ]\n";
			last;
		}

		# save any parsed data for this domain or barf loudly
		my $cnt_PIs_in_HMM = scalar keys %domain_data_per_pi;
		if (!$cnt_PIs_in_HMM)
		{
			push(@$HMMs_sans_hits, "$HMMoutput_file\t$HMM_name\t$HMM_acc");
#			print STDERR	"\nWarning!: No matching sequences in [ $HMMoutput_file ] "
#							."for [ $HMM_name : $HMM_acc ]\n";
		}
		else
		{
			# get the domain code for this HMM
			my @dom_code_thresholds = get_dom_code_thresholds($HMM_name, $HMM_acc, $HMM_filename, $dbh);
			if (@dom_code_thresholds)
			{
				save_domain_data_in_db($HMM_name, @dom_code_thresholds, %domain_data_per_pi, $dbh);
				$dom_codes{$dom_code_thresholds[0]}=1;
				print STDERR ":";
				print "$HMMoutput_file\t$HMM_name\t$HMM_acc\t$dom_code_thresholds[0]\n"	;
			}
			else
			{
				push(@$obsolete_hmms, $HMM_filename."\t".$HMM_name."\t[".$HMM_acc.']');
			}
			

		}

		# skip to end of record
		while (<HMMOUT>)
		{
			last if m{^//};
		}
	
	}

	close HMMOUT;

	return scalar keys %dom_codes;
}


#_________________________________________________________________________________________

#	get_dom_code_thresholds

#		get domain code corresponding to the hmm name, accession and filename

#_________________________________________________________________________________________
sub query_sql_for_single_tuple($$$\@)
{
	my ($dbh, $sql_str, $warn_of_multiple_tuples, $ref_results) = @_;
	my $sth	= $dbh->prepare($sql_str);
	$sth->execute() or die "\nError:\n$DBI::errstr\n";
	
	my $matches = $sth->rows();
	if ($matches > 1)
	{
		warn $warn_of_multiple_tuples;
		return $matches;
	}
	if ($matches == 1)
	{
		@$ref_results = $sth->fetchrow_array();
		$sth->finish();

		return $matches;
	}
	return 0;

}

sub get_dom_code_thresholds($$$$)
{
	my ($HMM_name, $HMM_acc, $HMM_filename, $dbh) = @_;

	my $err_msg = "Error!\nThe HMM Search output owner cannot be determined:\n"
					."\tMore than one domain has the name [ $HMM_name ], "
					."\taccession [ $HMM_acc ] "
					."\tand file name [ $HMM_filename ].\n";

	# get the domain code corresponding to this entry
	my @dom_code_thresholds;
	my $matches = query_sql_for_single_tuple($dbh,
									 'SELECT dom_code, thres_true, thres_maybe '
									.'FROM hmms natural join domains '
									.'WHERE '
										."hmm_name = '$HMM_name' AND "
										."hmm_acc = '$HMM_acc' ",
							$err_msg, @dom_code_thresholds);
	return @dom_code_thresholds unless ($matches == 0);

	# no matches for name and filename: try again
	$matches = query_sql_for_single_tuple($dbh,
									 'SELECT dom_code, thres_true, thres_maybe '
									.'FROM hmms natural join domains '
									.'WHERE '
										."hmm_file = '$HMM_filename' ",
							$err_msg, @dom_code_thresholds);

	return @dom_code_thresholds;
}	


#_________________________________________________________________________________________

#	save_domain_data_in_db($\@\%$)

#_________________________________________________________________________________________
sub save_domain_data_in_db($\@\%$)
{
	my ($HMM_name, $dom_code_thresholds, $domain_data_per_pi, $dbh) = @_;
	
									

	# delete all from domains_per_protein for this HMM
#	$dbh->do('DELETE FROM domains_per_protein WHERE dom_code = ?',
#				undef, $dom_code_thresholds->[0])  or die "\nError:\n$DBI::errstr\n";

	
	
	foreach my $pi (keys %$domain_data_per_pi)
	{
		my ($start, $finish, $p_value) = @{$domain_data_per_pi->{$pi}};
		my ($dom_code, $thres_true, $thres_maybe) = @$dom_code_thresholds;

		# Is this over the maybe or the true threshold?
		my $confidence;
		if ($p_value < $thres_true)
		{
			$confidence = 1;
		}
		elsif ($p_value < $thres_maybe)
		{
			$confidence = 0;
		}
		else
		{
			next;
		}

		# get two pi numbers
		my ($pi1, $pi2) =split /_/, $pi;
		my @data;
		push (@data,
						$pi1,
						$pi2,
						$dom_code,
						$HMM_name,
						$start,
						$finish,
						$p_value,
						$confidence);
		$dbh->do( 	'INSERT into domains_per_protein VALUES ('.
					join (	", ", ('?') x @data).
					')', undef, @data)
			or die "\nError:\nCould not insert data into domains_per_protein table\n".
					"pi=($pi1, $pi2), hmm=($HMM_name, $dom_code) ".
					(join ",", @{$domain_data_per_pi->{$pi}}).
					"\n$DBI::errstr\n";
		$cnt_dom_prot++;
	}
}

#_________________________________________________________________________________________

#	fasta_to_pi_callback($$$$)

#_________________________________________________________________________________________
sub fasta_to_pi_callback($$$$)
{
	# parameters
	my ($accession, $sequences, $line_num, $seq_cnt, $acc_to_pi) = @_;

	# print progress
	print STDERR '.' if (!($seq_cnt % 5000));

	unless ($accession =~ /([^|:]+[:|][^|]+)\|/)	#/([^|]+\|[^|]+\|)/)
	{
		print STDERR "\nError:\nAccession not found in $accession!!\n",
					 "at line=$line_num, record number = $seq_cnt\n";
		return;
	}

	# save pi and accession link
	my ($crc_lo, $crc_hi) = crc (join('', @$sequences), 64);
	$acc_to_pi->[($seq_cnt-1)*2] = $1;
	$acc_to_pi->[($seq_cnt-1)*2+1] = $crc_lo.'_'.$crc_hi;
}



#_________________________________________________________________________________________

#	create_postgres_tables


#_________________________________________________________________________________________
sub create_postgres_tables($)
{
	my ($dbh) = @_;	
	# this procedure tends to fail (if the tables already exist)
	# so no error messages but auto commit
	my @original_dbh_flags = @{$dbh}{'PrintError', 'AutoCommit', 'RaiseError'};
	@{$dbh}{'PrintError', 'AutoCommit', 'RaiseError'} = (0, 1, 0);

	if 	($dbh->do( 'CREATE TABLE domains_per_protein ('.
					'pi1               OID,'.
					'pi2               OID,'.
					'dom_code          INTEGER,'.
					'dom_name          TEXT,'.
					'start             INTEGER,'.
					'finish            INTEGER,'.
					'p_value           DOUBLE PRECISION,'.
					'confidence        INTEGER)'))
	{
		$dbh->do( 'CREATE INDEX i_domprot__pi on domains_per_protein(pi1,pi2)');
		$dbh->do( 'CREATE INDEX i_domprot__dom_code on domains_per_protein(dom_code)');
		$dbh->do( 'CREATE INDEX i_domprot__dom_name on domains_per_protein(dom_name)');
	}

	@{$dbh}{'PrintError', 'AutoCommit', 'RaiseError'} = @original_dbh_flags;
}	



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Core logic!!

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# connect to database
$| = 1;
my $dbh = DBI->connect( "dbi:Pg:dbname=$dbname",
                        $username,
                        $userpassword,
						{
                          RaiseError => 0,
                          AutoCommit => 0,
						  PrintError => 0
                        }
                      )
		or die "Fatal Error: Database connection could not be made: $DBI::errstr";

# make sure HMM / domain tables are present in the database
create_postgres_tables($dbh);

#drop all indices on the tables
my @sql_recreate_indices = drop_indices($dbh, 'domains_per_protein');


# save HMMoutput files
my @HMM_out_files = @ARGV;

# map accession to panda identifiers
my %acc_to_pi;
print STDERR "Parsing the fasta file:\n[";
@ARGV = ($seq_file);
my $cnt_seq;
eval
{
	my @acc_to_pi;
	$#acc_to_pi = 1_000_000;
	$cnt_seq = parse_fasta(\&fasta_to_pi_callback, \@acc_to_pi);
	$#acc_to_pi = $cnt_seq * 2 - 1;
	%acc_to_pi = @acc_to_pi;
};
print STDERR "]\n", "$cnt_seq sequences parsed.\n";


my %domain_cnts;
my @obsolete_hmms;
my @HMMs_sans_hits;
print STDERR "[";
foreach my $HMM_out_file (@HMM_out_files)
{
	my $domains_cnt_in_file =
		# calls iteratively save_domain_data_in_db($$\%$)
		parse_HMM_output_file($HMM_out_file, $seq_file, %acc_to_pi, $dbh,
							  $cnt_seq, @obsolete_hmms, @HMMs_sans_hits);
	$domain_cnts{$domains_cnt_in_file}++;
}
print STDERR "]\n";

#recreate indices
$dbh->do($_) or die "\nError:\n$DBI::errstr\n" for (@sql_recreate_indices);

$dbh->commit();
$dbh->disconnect();
foreach my $cnts(keys %domain_cnts)
{
	print STDERR $domain_cnts{$cnts},"\tfiles(s) contained matches from $cnts \thmm(s) each.\n" if $cnts;
}

if ($domain_cnts{0})
{
	print STDERR $domain_cnts{0},"\tfiles(s) did not contain any matches from a relevant hmm.\n";
}
# list of hmms which are obsolete
if (@obsolete_hmms)
{
	print STDERR scalar @obsolete_hmms, " HMM records do not have ".
			"corresponding data in the SMART database.\n\n";
	print "These domains PFAM are most probably obsolete and will not be entered into panda:\n";
	print join "\n", @obsolete_hmms, "", "";
}
if (@HMMs_sans_hits)
{
	print "These HMMs did not have hits for the specified sequences:\n";
	print join "\n", @HMMs_sans_hits, "", "";
}

print "\n\n$cnt_dom_prot domain predictions parsed and saved.\n"
STDOUT->autoflush(1);
_exit(0);

#./panda_enter_PFAM_HMM_into_postgres.pl -a /home/share/database/pfam/Pfam-A.seed -hmm /home/share/database/pfam/pfam.hmm > obsolete_pfam_hmms.res
#./panda_enter_SMART_HMM_into_postgres.pl -s=panda /home/share/database/smart/*.HMM > obsolete_smart_hmms.res
#./panda_parse_HMMOUT_into_db.pl  -s /home/share/database/nr/nr_smart /home/share/database/hmmout/*.hmmout

